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Introduction 


This  work  is  concerned  with  theoretical  methods  for  designing  individualized  optimal 
strategies  of  breast  cancer  surveillance.  The  problem  of  optimal  cancer  surveillance 
is  set  up  as  a  search  for  optimal  scheduling  of  screening  examinations  subject  to 
certain  constraints  on  the  number  and  timing  of  medical  tests.  The  hypothesis  to  be 
tested  is  that  the  efficacy  of  breast  cancer  detection  can  be  enhanced  through  incor¬ 
porating  aggregated  family  history  information  into  a  mathematical  model  designed 
to  construct  optimal  schedules  of  cancer  surveillance.  The  proposed  methods  are  to 
be  validated  using  epidemiological  data  on  breast  cancer  from  the  Utah  Population 
Data  Base  and  the  Utah  Cancer  Registry. 


1.  Statement  of  Work 

This  annual  report  covers  the  following  four  tasks  formulated  in  the  statement  of 
work. 

Task  1:  Develop  programs  for  estimating  the  hazard  function  for  breast  cancer  using 
the  data-adaptive  method  of  kernel  estimation. 

Task  2:  Develop  programs  for  estimating  the  hazard  function  for  breast  cancer  using 
the  spline  approximation  method. 

Task  3:  Develop  programs  for  estimating  the  hazard  function  for  breast  cancer  using 
parametric  estimation  techniques. 

Task  4:  Construct  optimal  surveillance  scheduling  strategies  based  on  the  three 
methods  of  hazard  function  estimation. 


2.  The  research  carried  out  to  meet  the  objectives 
of  Tasks  1  and  2 

2.1.  Introduction 

The  shape  of  the  hazard  function  may  lead  to  insights  into  the  biology  of  carcino¬ 
genesis  which  may  not  be  easily  discernable  from  a  study  of  the  survival  function 
alone.  For  example,  it  is  typical  in  the  analysis  of  tumor  recurence  data  to  find  a 
hazard  function  that  is  bimodal  or  unimodal,  and  that  tends  to  zero  as  time  tends  to 
infinity  [1].  The  modes  of  the  hazard  may  be  interpreted  biologically  as  arising  from 
two  different  types  of  failure,  one  that  tends  to  occur  earlier  and  one  that  tends  to 
occur  later.  In  the  context  of  the  age-specific  hazard  function  for  cancer  incidence 
the  decrease  in  the  hazard  function  to  zero  may  lead  one  to  conclude  that  there  is 
a  non-zero  fraction  of  unaffected  (immune)  individuals.  In  fact,  if  the  cumulative 
hazard  appears  to  be  bounded,  one  should  expect  the  existence  of  a  non-zero  immune 
fraction.  More  generally,  a  large  degree  of  heterogeneity  in  disease  susceptibility 
may  lead  to  a  population  hazard  function  with  one  or  more  well-defined  maxima. 
The  maxima  may  correspond  to  discrete  subpopulations  with  different  genetic  pre¬ 
disposition  to  disease.  A  maximum  may  also  result  from  a  continuous  frailty,  as 
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the  surviving  population  at  higher  ages  may  be  overrepresented  by  individuals  with 
lower  risk  [2],  On  the  other  hand,  a  monotonically  increasing  (with  age  at  detection) 
hazard  rate  is  consistent  with  a  popular  belief  that  spontaneous  tumorigenesis  can 
be  modeled  as  a  Poisson  process. 

The  form  of  the  hazard  function  for  breast  cancer  incidence  may  depend  on 
characteristics  inherited  susceptibility.  Inherited  mutations  in  p53,  BRCA1,  BRCA2, 
the  ataxia-telangiectasia  gene  (AT),  HRAS,  and  the  androgen  receptor  gene  (AR) 
have  been  shown  to  play  a  role  in  breast  cancer  susceptibility  [3].  About  56%  of 
carriers  of  the  mutation  BRCA1  or  BRCA2  will  get  breast  cancer  by  the  age  of  70 
years  [4] .  BRCA1  has  an  estimated  allele  frequency  of  between  0.0002  and  0.001  (95% 
Cl)  [5],  and  accounts  for  about  3%  of  diagnosed  breast  cancer  [6].  The  allele  frequency 
of  mutations  in  BRCA2  is  estimated  at  0.00022  [7].  Germline  mutations  in  p53 
and  AR  are  extremely  rare,  and  mutations  in  the  HRAS1  minisatellite  locus  which 
confer  increased  risk  of  breast  cancer  are  also  rare,  having  an  estimated  population 
frequency  of  6%  [3].  In  a  study  of  100  Finnish  breast  cancer  families  analyzed 
by  protein  truncation  tests  and  direct  sequencing,  Vehmanen  et  al.  [8]  found  that 
only  21%  of  breast  cancer  families  were  accounted  for  by  mutations  of  BRCA1  and 
BRCA2,  providing  indirect  evidence  for  the  existence  of  other,  undiscovered  breast 
cancer  genes. 

Additional  insight  can  be  gleaned  from  the  hazard  function  for  cancer  incidence  in 
the  framework  of  a  mechanistic  model  of  carcinogenesis.  The  most  widely  accepted 
model  is  the  Moolgavkar-Venzon-Knudson  two-stage  clonal  expansion  model  [9,10]. 
The  Moolgavkar-Venzon-Knudson  model  has  the  following  assumptions: 

(A)  Normal,  susceptible  target  cells  are  initiated  according  to  a  (nonhomogeneous) 
Poisson  process  with  intensity  u(t). 

(B)  The  expansion  of  the  colony  of  initiated  cells  and  malignant  transformation  is 
specified  by  a  stochastic  birth-death-migration  process  with  the  division,  death  (or 
differentiation)  and  transformation.  Premalignant  cells  either  divide  into  two  pre- 
malignant  cells  with  rate  a(t),  die  with  rate  /3(f),  or  divide  asymmetrically  into  one 
premalignant  cell  and  one  malignant  cell  with  rate  p(f). 


It  has  been  shown  that  the  hazard  function  for  the  Moolgavkar-Venzon-Knudson 
model  with  constant  parameters  increases  monotonically  and  approaches  an  asymp¬ 
tote  [11].  An  asymptotic  value  for  the  hazard  is  also  reached  for  the  Moolgavkar- 
Venzon-Knudson  model  with  piecewise  constant  parameters,  and  in  that  case  the 
value  of  the  asymptote  depends  only  on  the  value  of  the  coefficients  in  the  unbounded 
interval  [11,12]. 

Expressions  for  the  survivor  function  were  first  obtained  by  Moolgavkar  and  Lue- 
beck  [11] .  A  simple  explicit  formula  for  the  survivor  function  S (t)  for  the  Moolgavkar- 
Venzon-Knudson  model  with  constant  parameters  was  obtained  by  Kopp-Schneider 
et  al.  [13]  and  Zheng  [14]: 


S(t ) 


2ce°'5(“a+/3+At“c)t 

(-a  +  [3  +  ij,  +  c)  +  (oi-(3-ij,  +  c)e~ct 


(1) 


where  c  =  <J(a  +  /3  +  ft)2  -  4a/3.  Zheng  also  presented  an  expression  for  the  proba- 
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bility  generating  function  for  the  number  of  malignant  cells  given  a  single  malignant 
cell  at  time  t  =  0,  allowing  an  expression  for  the  promotion  time  distribution 

_  (a  —  (d  —  jj,  +  c)(a  —  f3  —  (x  —  e)e~ct  +  {a  —  (3  —  //  +  c)(—a  +  P  +  y-  +  c) 
2a[(a  —  P  —  y  +  c)e~ct  +  (—a  +  ft  +  y  +  c)] 

(2) 

to  be  given.  It  is  easy  to  see  that  S(t)  and  F(t)  above  are  related  by  the  formula 

S(t)  =  exp  j— vj  F(x)dx  j  (3) 

which  was  shown  by  Hanin  and  Yakovlev  [15]  to  be  valid  in  a  more  general  setting. 

Yakovlev  and  Tsodikov  [16]  replace  assumption  (B)  above  with  the  following  as¬ 
sumption: 

(C)  Progenitor  cells  are  transformed  into  malignant  lesions  at  a  random  with  cumu¬ 
lative  distribution  function  F(x).  All  progenitor  cells  are  promoted  independently 
of  one  another. 

Assuming  F(0)  =  0,  it  follows  that  the  process  of  malignant  transformation  is 
also  a  Poisson  process,  with  integral  rate  A (t)  =  Jq  u{u)F{t  —  u)du.  As  in  the 
Moolgavkar-Venzon-Knudson  model,  the  simplest  model  of  spontaneous  carcinogenis 
takes  v(t)  =  v  to  be  constant,  in  which  case  A (t)  =  u /04  F{u)du  and  the  hazard 
function  for  time-to-tumor,  given  by  A (t)  =  vF{t),  is  nondecreasing.  The  probability 
S(t)  that  there  are  no  malignancies  by  time  t  is  then  given  by  (3). 

This  model  may  easily  be  modified  to  handle  inherited  lesions,  via  the  limiting 
case  where  v  is  taken  to  be  a  delta  function  at  the  origin.  If  F(t)  is  assumed  to  be 
absolutely  continuous,  then  the  integral  hazard  rate  A (t)  is  equal  to  uF(t )  and  the 
hazard  function  A (t)  =  uF'(t)  =  i where  f(t)  is  the  probability  density  function 
associated  with  F(t).  We  see  that  the  hazard  function  for  spontaneous  and  inherited 
lesions  are  quite  likely  to  have  very  different  shapes. 

Even  though  a  thorough  study  of  the  hazard  function  may  lead  to  new  insight 
into  the  process  of  carcinogenesis,  few  if  any  population-based  cohorts  have  been 
analyzed  to  determine  the  hazard  function  for  cancer  incidence. 

2.2.  Data 

The  data  for  this  study  was  obtained  by  linking  records  from  the  Utah  Population 
Data  Base  (UPDB)  with  the  Utah  Cancer  Registry  (UCR).  The  UPDB  consists  of 
the  genealogical  records  of  more  than  1,000,000  individuals  who  were  born,  died,  or 
married  in  Utah,  or  en  route  to  Utah  during  the  nineteenth  and  twentieth  centuries. 
Since  1973  the  UCR  has  been  reporting  to  National  Cancer  Institutes  Surveillance 
Epidemiology  and  End  Results  (SEER)  program,  and  is  required  to  maintain  very 
high  standards  for  case  reporting  and  follow-up,  and  to  periodically  undergo  quality 
control  audits  by  SEER  personnel  to  assure  uniformly  high  quality  and  consistency 
from  year  to  year.  The  available  follow-up  information  comes  either  from  Utah  death 
certificates,  which  have  been  linked  to  the  UPDB  genealogical  data  every  year  from 
1933  through  the  beginning  of  1997,  or  from  linkage  of  the  HCFA  beneficiary  data 
to  the  UPDB.  The  study  population  consists  of  126,141  men  and  122,208  women 
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recorded  in  the  Utah  Population  Database,  who  were  born  from  1874  to  1931  and  for 
whom  follow-up  information  is  available  that  places  them  in  Utah  during  the  years 
of  operation  of  the  Utah  Cancer  Registry  (1966-present).  There  are  5,372  cases  of 
female  breast  cancer  represented  in  the  data.  Since  methodological  research  was  the 
major  undertaking  in  Year  1,  we  plan  to  conduct  analyses  of  epidemiological  data 
in  2000.  Analyses  will  be  performed  on  subcohorts  based  on  birth  year  (1874-1889, 
1890-1899,  1900-1909,  1910-1919,  and  1920-1931). 


2.3.  Methodological  Problems 

Truncation:  Nonparametric  Estimation 

We  wish  to  estimate  the  age  specific  hazard  function  for  breast  cancer  from  the  data 
described  above,  taking  into  account  that  the  data  is  subject  to  random  truncation: 
cases  which  occurred  during  or  before  1965  are  not  recorded  in  the  dataset.  Subject 
were  between  the  ages  of  34  and  86,  at  the  time  of  truncation.  Thus,  analysis  of  the 
data  must  take  into  account  not  only  to  the  effects  of  right  censoring,  but  also  the 
effects  of  left  random  truncation  due  to  delayed  entry  into  the  risk  set. 

The  problem  of  random  truncation  can  be  formulated  as  follows.  Let  the  trun¬ 
cation  time  Y  have  distribution  function  G(y)  and  the  failure  time  (time  of  cancer 
diagnosis)  X  have  distribution  function  F(x).  We  require  that  truncation  be  indepen¬ 
dent  of  failure  and  for  simplicity  assume  no  censoring  for  the  present.  Observations 
are  conditional  on  the  random  event  X  >  Y.  Let  G*(y )  and  F*{x)  be  the  corre¬ 
sponding  distribution  functions,  conditional  on  X  >  Y.  Let  S(x)  =  1  —  F(x)  be  the 
survivor  function  of  X.  Suppose  that  we  have  observations  (Yj*, X±), . . . ,  (Y*,X*) 
from  the  conditional  distribution.  The  full  likelihood  of  the  observed  data  is  given 

by 

L=  n  [dFiX^dGW/a] ,  (4) 

i= i 

where  a  =  f  fy<x  dF(x)dG(y).  A  key  observation  is  that  if  X  and  Y  are  independent, 
then  the  hazard  of  X  given  X  >  Y  =  y  a,t  x  >  y  is  equal  to  the  hazard  of  X  at  x 
[17,18].  This  observations  leads  to  the  result,  first  mentioned  by  Kaplan  and  Meier 
[19],  that  if  the  distribution  G{t)  is  allowed  to  vary  freely,  the  natural  generalization 
of  the  product  limit  estimator,  given  by  the  formula 


S(t)  = 


(5) 


where  R(u)  =  #{Y*  <  U  <  X*}  is  the  number  of  subjects  at  risk  at  U,  is  the 
nonparametric  maximum  likelihood  estimator  (NPMLE)  of  the  survivor  function 
S(t)  of  X  (see,  for  example  [17,18,20]). 

This  result  extends  naturally  to  the  case  with  random  independent  censoring  [20]. 
It  also  easily  follows  that  in  the  nonparametric  setting  (again  with  no  censoring), 
maximizing  (1)  is  equivalent  to  maximizing  the  conditional  likelihood  of  (Xj* , . . . ,  X*) 
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given  (Tj*, . ..  ,Y*),  which  can  be  written 


CL  =  I]  f(X’)/S(Y').  (6) 

i=  1 

(see,  for  example,  [19-22]).  Maximizing  the  conditional  likelihood  also  leads  to  the 
familiar  Nelson-Aalen  estimator  for  the  integrated  hazard  function  H(t)  of  X  [20], 
which  is  given  by 

A(f)  =  £  R(X>)~\  (7) 

X*<t 

These  results  can  be  extended  to  the  case  of  right  censoring  [20] . 


Truncation:  Parametric  Models 

We  consider  the  situation  where  X  and  Y  are  independent,  F(x)  is  parametrized, 
while  G(y )  is  allowed  to  vary  freely.  In  a  later  subsection  F(x)  will  be  come  from  a 
quadratic  spline  model. 

The  data  are  independent  pairs  (yi,aq), . . . ,  (yn,  xn)  from  the  joint  distribution 
(Y,X),  conditional  on  (Y  <  X).  We  suppose,  for  simplicity,  that  there  are  no  ties 
among  yi,  y2,  •  •  •  ,  Vn,  and  suppose  X  has  absolutely  continuous  distribution  function 
coming  from  a  family  F(x ;  z)  parameterized  by  a  vector  z,  with  corresponding  sur¬ 
vival  function  S(x',  z)  =  1  -  F(x ;  z)  and  density  f(x]  z).  The  NPMLE  for  G  should 
consist  of  (unknown)  point  masses  <?i,  <?2,  •  •  •  >  Qn  placed  at  the  points  yi,y2,  ■  ■  ■  ,Vn- 
The  logarithm  of  the  complete  likelihood  (4)  can  be  rewritten 


log(L)  =  ^[log(/(a:i;  z))  +  log(^)]  -  nlog 

i—1 


J2S(yl-,z)qi 


j=l 


(8) 


If  we  factor  the  out  the  part  of  the  likelihood  corresponding  to  formula  (6),  the 
logarithm  is  given  by 


log(CL)  =  ?))  ~  l°g(s(,yu  z))]-  (9) 

i=l 

We  now  discuss  the  changes  which  must  be  made  in  when  censoring  and  additional 
covariates  are  present.  If  s  is  a  vector  of  additional  covariates,  A  (a:,  s;  z)  denotes  the 
hazard  associated  with  F(x,  s;  z)  and  A(x,  s;  z)  the  cumulative  hazard,  we  note  that 
(9)  becomes 


log(CL)  =  ^)[log(A(a:i,  s;  z))  -  (A(xi,  s);  z)  -  A  (yi,  s;  z))].  (10) 

i= 1 

In  the  presence  of  right  censoring  which  is  independent  of  both  the  failure  and 
truncation  times,  aq  is  replaced  in  the  above  formulation  by  the  minimum  of  the 
failure  and  censoring  time.  The  term  f{x,s\z )  in  the  likelihood  is  replaced  by 
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f(x,  s ;  z)sS(x ,  s;  z)1  s,  where  <5*  =  1  if  observation  i  is  a  failure  and  A  =  0  otherwise, 
and  the  conditional  likelihood  (6)  (with  Xi,Si  and  yi  regarded  as  fixed)  becomes 

CL  =  f\[f{xhSi\Zi)& S{xuSi\Zi){1~6)}/ S{yi,Si\Zi). 

2=1 

In  this  setting  log(C'L)  becomes 

n 

loglC*-!!)  =  )  ',\Sj  k)§(A(Xj,  ^i))  (A(Xj,  Sj,  z'j  A(yi,  Si ,  2^))]  ■  (11-) 

i= 1 

W  havee  chosen  to  maximize  (11)  rather  than  the  full  likelihood.  Under  appropriate 
conditions  the  two  maximization  procedures  are  equivalent. 

Spline  Models 

We  choose  to  model  the  hazard  via  quadratic  splines  as  in  [23].  A  quadratic  spline 
with  m  knots  specifies  the  hazard  to  be  of  the  form 

2  m 

Xm(t)  =  5Z  70 +  J2  7?2(*  -  U)+  (12) 

2=0  j=l 

where  (x)+  =  max(x,0).  For  each  birth  cohort,  we  will  fit  splines  with  knots  which 
are  equally  spaced  in  the  interior  of  the  interior  [Tmin,Tmax],  where  Tmin  is  the  min¬ 
imum  truncation  age  in  the  cohort  and  Tmax  the  maximum  follow-up  (failure  or 
censoring)  time.  Restrictions  should  be  placed  on  the  coefficients  to  ensure  that  A m(t) 
remains  positive  for  all  t.  Thus  with  m  knots  the  number  of  parameters  is  m  +  3. 
Models  can  be  fit  using  maximum  likelihood  techniques  applied  to  the  conditional 
likelihood,  as  given  by  (11). 

We  have  developed  software  designed  to  compute  the  spline  estimates  by  max¬ 
imizing  log(CL)  using  the  algorithm  of  Powell  [24].  We  start  with  one  knot  and 
increase  the  number  of  knots  until  the  fit  is  not  improved,  as  determined  by  the  like¬ 
lihood  ratio  test  at  the  significance  level  a  =  0.05.  Two  other  subcohort  estimates 
of  the  hazard  function  can  be  computed  for  comparison  with  the  spline  estimator;  a 
life  table  version  of  (5),  and  a  Gaussian  kernel  estimate  based  on  the  Nelson- Aalen 
estimator  (7).  We  have  developed  relevant  software  for  this  purpose  as  well. 

2.4.  Future  Plans 

Using  the  computer  programs  developed  in  Year  1,  the  hazard  function  for  cancer 
incidence  will  be  estimated  from  left  truncated  and  right  censored  data  based  on 
the  the  conditional  likelihood.  Four  estimation  procedures  based  on  the  conditional 
likelihood  will  be  used  to  estimate  the  age-specific  hazard  function  from  the  data; 
these  are  the  life-table  method,  a  kernel  method  based  on  the  Nelson  Aalen  estimator, 
a  spline  estimate,  and  a  proportional  hazards  estimate  based  on  splines  with  birth 
year  as  sole  covariate.  The  latter  model  is  aimed  at  adjusting  for  the  cohort  effect. 
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3.  The  research  carried  out  to  meet  the  objectives 
of  Tasks  3  and  4 


3.1.  Introduction 

When  evaluating  potential  benefits  of  screening  in  different  populations  (different 
racial  composition,  different  geographic  regions,  etc.),  it  makes  sense  to  compare  the 
expected  effects  for  the  schedules  which  are,  in  some  clearly  defined  sense,  optimal  for 
each  of  the  populations  under  consideration;  this  represents  the  usual  way  of  "stan¬ 
dardization  through  optimization1 2'.  In  doing  so,  we  compare  maxima  of  potential 
benefits  that  could  be  gained  in  each  setting. 

The  problem  of  screening  optimization  is  of  great  interest  in  itself.  Because  of 
the  significant  cancer  incidence  and  progress  of  tumor  detection  technology,  cancer 
surveillance  and  screening  are  becoming  increasingly  important  and  costly  public 
health  problems.  It  is  clear  that  appropriate  mathematical  methods  are  indispens¬ 
able  for  a  more  effective  management  of  the  caseload  through  designing  optimal 
surveillance  strategies.  Interest  in  exploring  this  avenue  has  quickened  in  the  past 
few  decades  [16,  25-39]. 

The  problem  of  optimal  cancer  surveillance  can  be  set  up  as  a  search  for  opti¬ 
mal  scheduling  of  screens  subject  to  certain  constraints  on  the  number  and  timing 
of  medical  examinations.  Yakovlev  and  Tsodikov  [16]  have  developed  methods  for 
constructing  optimal  surveillance  strategies  based  on  the  minimum  delay  time  crite¬ 
rion,  given  that  the  total  number  of  examinations  is  fixed,  see  also  [39].  They  used 
dynamic  programming  methodology  to  solve  the  associated  optimization  problem. 
However,  there  are  two  weak  points  in  their  approach.  First,  the  probability  of  tumor 
detection  is  assumed  to  be  independent  of  the  process  of  tumor  regrowth.  Second, 
estimation  of  the  tumor  onset  time  distribution  is  feasible  only  if  a  sample  of  diag¬ 
nostic  times  produced  by  a  discrete  surveillance  program  with  known  false  negative 
rate  is  available.  The  same  applies  equally  to  pre-diagnosis  screening  programs. 

To  surmount  the  above  mentioned  problems  we  have  explored  a  new  avenue  in 
the  problem  of  cancer  screening  optimization.  The  newly  developed  methodology  is 
well  adapted  to  the  structure  of  the  data  amassed  in  the  Utah  Population  Data  Base 
and  the  Utah  Cancer  Registry.  As  a  measure  of  the  effect  of  screening,  we  propose  to 
use  the  difference  between  the  expected  tumor  sizes  at  detection  with  and  without 
screening,  which  coincides  with  the  Kantorovich  distance  between  the  distributions 
of  the  corresponding  random  variables.  The  structure  of  this  distance  allows  for 
characterizing  the  net  effect  of  screening,  as  compared  to  that  of  spontaneous  detec¬ 
tion.  Taken  alone,  the  design  of  optimal  schedules  of  cancer  surveillance  does  not 
require  information  on  tumor  size  at  detection;  such  schedules  can  be  constructed 
once  the  basic  model  parameters  have  been  estimated  in  one  way  or  another  from 
epidemiological  data  at  hand. 

The  proposed  approach  offers  the  following  distinct  advantages: 

1.  It  provides  a  simple  but  still  realistic  description  of  cancer  latency; 

2.  It  can  be  generalized  in  various  ways  while  retaining  its  basic  structure; 
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3.  It  furnishes  a  biologically  meaningful  interpretation  of  data  analyses; 

4.  It  accommodates  standard  population-based  statistical  data;  its  implementation 
does  not  depend  heavily  on  availability  of  the  data  yielded  by  screening  trials; 

5.  Rigorous  statistical  methods  are  available  for  estimation  of  model  parameters; 

6.  It  can  be  used  for  designing  optimal  strategies  of  cancer  screening  and  surveillance. 

3.2.  The  Model 

In  describing  the  natural  history  of  cancer,  the  process  of  tumor  development  can  be 
broken  down  into  three  stages.  These  stages  are: 

•  formation  of  initiated  cells; 

•  promotion  of  initiated  cells  resulting  in  appearance  of  the  first  malignant  clono- 
genic  cell; 

•  subsequent  growth  and  progression  of  malignant  tumor. 

The  duration  of  each  stage  of  carcinogenesis  is  thought  of  as  a  random  variable. 
In  our  sample  calculations,  we  used  a  two-parameter  gamma  family  to  specify  the 
distribution  of  the  length  of  the  first  two  stages  of  carcinogenesis.  However,  more 
elaborate  mechanistic  models  of  carcinogenesis  are  available  to  describe  the  time 
to  the  event  of  malignant  transformation  at  the  cellular  level  (see  Section  2.1).  In 
particular,  we  plan  to  use  the  Moolgavkar-Venzon-Knudson  model  and  the  Yakovlev- 
Polig  model  in  our  future  research. 

We  proceed  from  the  following  general  functional  form  for  the  tumor  size  (the 
number  of  cells  in  a  tumor)  S  : 


S(w)  =  fe(w), 

where  w  is  the  time  from  the  moment  of  the  onset  of  cancer,  and  9  is  a  parameter 
which  may  be  scalar  or  vector,  deterministic  or  random.  It  is  assumed  that,  for  every 
9,  fe  is  a  strictly  monotonically  increasing  absolutely  continuous  function  such  that 
fd( 0)  =  1.  For  a  given  9,  denote  by  ge  the  inverse  function  for  fe,  and  set 

rw 

$o(w)  :=  /  fe(u)du. 

Jo 

Specific  laws  of  tumor  growth  of  primary  interest  are: 

(1)  Deterministic  exponential  growth;  in  this  case,  S(w )  =  e*w,  where  A  >  0  is  a 
constant  growth  rate,  see  [40]  for  substantiation; 

(2)  Exponential  growth  with  A  thought  of  as  a  gamma  distributed  r.  v.  [41]; 

(3)  The  Gompertz  law: 

S(w)  =  eA^e-Bw\ 
with  constant  parameters  A,  B  >  0. 
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The  sequence  of  moments  of  time  assigned  for  medical  exams  for  a  specific  cancer 
and  counted  from  the  birth  of  a  patient  will  be  called  a  screening  schedule.  Let  T  be 
the  set  of  all  possible  screening  schedules  r  =  {t\  <  <  ...  <  rn}.  The  set  T  may 

be  subject  to  (some  of)  the  following  restrictions: 

(a)  n  <  no,  where  n0  is  an  upper  bound  for  the  number  of  exams; 

(b)  Tx>m  and  rn  <  M,  where  m  and  M  are  the  earliest  and  the  latest  times  for 
the  first  and  the  last  exams,  respectively; 

(c)  Ti+i  -  Ti  >  h  >  0  for  alH  =  1,2,  ...,n  -  1.  This  condition  suggests  a  lower 
bound  h  for  the  minimal  duration  between  any  two  successive  exams. 

Other  restrictions  on  the  moments  of  exams  can  also  be  accommodated.  In  the 
language  of  control  theory,  the  set  T  is  referred  to  as  the  set  of  admissible  schedules. 

3.3.  The  Screening  Efficiency  Functional 

Numerous  attempts  have  been  made  to  relate  the  probability  of  detecting  a  tumor 
to  its  size  [40-46].  Following  Brown  et  al.  [44],  we  assume  that  the  rate  of  tumor 
detection  is  proportional  to  the  current  tumor  size. 

We  distinguish  between  spontaneous  and  screening  based  tumor  detections.  The 
first  occurs  in  the  absence  of  or  concurrently  with  screening  and  is  thought  of  as  a 
continuous  process.  In  contrast  to  this,  screening  based  detection  is  an  instantaneous 
event  that  may  occur  only  at  the  moments  of  the  prescribed  medical  exams  and  is 
therefore  a  discrete  process.  When  both  types  of  detection  are  present,  they  can  be 
viewed  as  competing  risks. 

Let  random  variables  W0  and  W1}  denote  the  times  of  spontaneous  and  screening 
based  detections,  measured  from  the  moment  of  cancer  onset,  respectively,  and  let 
T  be  the  time  of  tumor  onset.  We  have  derived  a  formula  for  the  screening  efficiency 
functional  proceeding  from  the  following  two  biologically  natural  assumptions. 

1.  The  r.v.’s  W0  and  T  are  independent. 

2.  For  every  t  >  0,  the  r.v.’s  Wa  and  W0  are  conditionally  independent  given  that 
T  =  t. 

The  first  assumption  implies  that  the  moment  of  spontaneous  tumor  detection  mea¬ 
sured  from  the  appearance  of  the  first  malignant  clonogenic  cell  is  independent  of 
the  prior  duration  of  tumor  latency.  The  second  assumption  reflects  a  technological 
(or  instrumental)  nature  of  both  detection  processes.  It  states  that,  given  the  mo¬ 
ment  of  cancer  onset,  the  two  times  W0  and  Wlt  at  which  competing  events  of  the 
spontaneous  and  screening  based  tumor  detection  may  occur,  are  independent.  This 
statement  immediately  follows  from  the  assumption  that  both  detection  processes 
are  completely  determined  by  the  current  tumor  size  as  a  deterministic  function  of 
time. 

For  an  admissible  screening  schedule  r  €  T,  we  define  the  efficiency  functional 
as  the  Kantorovich  distance  dK(N0,N;T )  between  the  tumor  sizes  N0  and  N  at 
spontaneous  and  combined  detection.  It  is  well  known  [47,  48]  that 

/OC  _ 

|  FNo(n)  -  FN(n)  |  dn. 
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An  alternative  expression  for  the  efficiency  functional  is  given  by 

/oo  _  roc  _ 

FN0(n)dn  —  J  Fif{n)dn  —  EN0  —  EN, 
where  E  stands  for  the  expectation. 

An  explicit  analytic  expression  for  screening  efficiency  functional  is  presented  in 
the  attached  paper  by  Hanin  et  al.  (accepted  for  publication  in  Mathematical  and 
Computer  Modelling). 

The  optimization  problem 


d(N,  AT0;  t )  —>  max,  tgT, 

has  been  solved  by  exhaustive  search  with  some  simplification  arising  from  the  special 
form  of  the  dependence  of  the  efficiency  functional  on  r. 

3.4.  Numerical  Experiments 

The  purpose  of  our  numerical  experiments  was  to  check  feasibility  of  numerical  and 
optimization  problems  associated  with  the  proposed  approach  to  stochastic  modeling 
of  cancer  screening.  These  experiments  are  described  at  length  in  the  attached  paper 
by  Hanin  et  al.  The  most  interesting  finding  was  that  in  the  case  of  exponential 
tumor  growth  optimal  screening  schedules  are  uniform  or  very  close  to  such.  The 
parameter  values  used  in  our  sample  computations  were  judiciously  chosen  on  the 
basis  of  our  experience;  they  are  no  better  than  an  educated  guess.  However,  this 
study  clearly  shows  that  mathematical  and  computational  problems  of  optimal  cancer 
surveillance  are  tractable  within  the  framework  of  the  proposed  model.  What  is  now 
required  is  to  make  a  final  stride  towards  the  use  of  parameter  estimates  obtained 
from  population-based  epidemiologic  data. 

3.5.  Future  Plans 

The  model  will  be  validated  using  the  relevant  epidemiological  data  on  birth  co¬ 
horts  identified  through  the  Utah  Population  Data  Base.  An  appealing  possibility 
would  be  to  estimate  unknown  parameters  of  breast  cancer  latency  solely  from  the 
age-at-diagnosis  data.  However,  this  is  unlikely  to  be  feasible  given  the  number  of 
parameters  incorporated  into  the  model.  The  main  obstacle  is  that  the  available  in¬ 
formation  is  too  sparse  for  estimating  all  the  parameters,  and  a  search  for  additional 
sources  of  data  is  clearly  warranted.  On  the  other  hand,  our  modeling  techniques 
allow  incorporation  of  the  process  of  tumor  detection  into  mechanistic  models  of 
tumor  latency,  thereby  making  it  possible  to  utilize  data  on  tumor  size  at  detec¬ 
tion  as  an  additional  source  of  information  on  the  natural  history  of  the  disease. 
This  information  is  available  in  the  Utah  Cancer  Registry.  For  a  given  functional 
form  of  tumor  growth,  we  will  obtain  within  the  framework  of  our  model  the  cor¬ 
responding  joint  distribution  of  tumor  size  and  age  at  detection.  Proceeding  from 
this  joint  distribution  we  will  construct  the  likelihood  of  the  sample  under  study. 
The  likelihood  function  will  be  maximized  numerically  by  computer  using  nonlinear 
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programming  methods  [24],  Once  the  estimation  problem  has  been  solved  we  are 
in  a  position  to  design  optimal  schedules  of  breast  cancer  surveillance  allowing  for 
individual  information  on  family  history. 


Key  Research  Accomplishments 

Our  key  accomplishments  in  Year  1  can  be  summarized  briefly  as  follows: 

•  We  have  developed  computer  programs  implementing  four  statistical  procedures 
for  estimation  of  the  hazard  function;  these  procedures  accommodate  data  subjected 
to  random  truncation  and  censoring. 

•  A  new  method  has  been  developed  for  designing  optimal  schedules  of  breast 
cancer  surveillance  specially  adapted  to  population-based  settings. 

•  Numerical  experiments  have  shown  that  mathematical  and  computational  prob¬ 
lems  of  optimal  cancer  surveillance  are  tractable  within  the  framework  of  the  pro¬ 
posed  model  of  cancer  surveillance  and  detection. 


Reportable  Outcomes 

1.  Hanin,  L.G.,  Tsodikov,  A.D.,  and  Yakovlev,  A.Y.  Optimal  schedules  of  cancer 
surveillance  and  tumor  size  at  detection,  Mathematical  and  Computer  Modelling ,  in 
press  (see  Appendix). 

2.  Hanin,  L.G.  Optimal  schedules  of  cancer  surveillance,  presented  at  the  Interna¬ 
tional  Workshop  on  Cure  Rate  Estimation,  Tampa,  Florida,  February  19-21,  1999. 

3.  Yakovlev,  A.Y.,  Tsodikov,  A.D.,  and  Hanin,  L.G.  Optimal  schedules  of  breast 
cancer  surveillance,  Abstract,  Era  of  Hope  Meeting,  Atlanta,  June  2000. 

4.  Boucher,  K.M.  and  Kerber,  R.A.  The  shape  of  the  hazard  function  for  cancer 
incidence,  Abstract,  Era  of  Hope  Meeting,  Atlanta,  June  2000. 

5.  Yakovlev,  A.Y.  Mechanistic  Modeling  of  Breast  Cancer  Surveillance,  Grant  Appli¬ 
cation,  RFA  "Cancer  Intervention  and  Surveillance  Network  (CISNET)",  NIH/NCI. 


Conclusions 

We  have  developed  a  mathematical  model  yielding  an  algorithm  for  designing  opti¬ 
mal  schedules  of  breast  cancer  surveillance.  An  explicit  expression  of  the  screening 
efficiency  functional  has  been  derived.  The  derivation  is  based  on  a  plausible  assump¬ 
tion  that  the  intensity  of  detection  (the  hazard  function  for  the  age  at  detection) 
is  proportional  to  the  current  tumor  size.  The  main  advantage  of  the  proposed  ap¬ 
proach  is  that  it  accommodates  cohort  data  of  a  fairly  general  structure,  not  only  the 
data  resulting  from  screening  trials.  We  have  also  developed  several  numerical  algo¬ 
rithms  and  software  for  estimating  the  hazard  function  for  breast  cancer  incidence 
from  the  data  amassed  in  the  Utah  Population  Data  Base.  Allowing  for  the  effects 
of  random  censoring  and  truncation,  these  procedures  will  be  used  in  future  research 


15 


for  testing  covariate  effects  associated  with  different  indicators  of  family  history. 


So  what?  We  now  have  the  necessary  tools  for: 

•  designing  optimal  schedules  of  breast  cancer  surveillance  given  the  numerical 
parameters  describing  the  natural  history  of  the  disease  are  known; 

•  testing  covariate  effects  associated  with  different  indicators  designed  to  aggre¬ 
gate  family  history  information. 

In  Year  2,  our  focus  will  be  on  the  development  of  methods  for  estimation  of  the 
natural  history  of  breast  cancer  from  the  data  available  in  the  Utah  Population  Data 
Base  and  the  Utah  Cancer  Registry.  Using  these  resources  we  will  also  test  indicators 
of  family  history  in  order  to  select  the  most  informative  one  for  the  purposes  of 
individualization  of  breast  cancer  surveillance  strategies. 
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ABSTRACT  -  The  paper  explores  methodological  and  mathematical  aspects  of  a 
new  approach  to  constructing  optimal  schedules  of  cancer  screening.  This  approach 
consists  in  systematic  use  of  tumor  size  at  detection,  combining  stochastic  models  of 
tumor  latency,  tumor  growth  and  tumor  detection,  and  employing  a  new  biologically 
natural  screening  efficiency  criterion  defined  as  the  Kantorovich  distance  between  the 
tumor  size  at  spontaneous  detection  in  the  absence  of  screening  and  the  tumor  size  at 
detection  when  both  spontaneous  and  screening  based  mechanisms  are  in  place.  An 
explicit  formula  for  the  efficiency  functional  is  obtained.  Sample  calculations  suggest 
that  in  the  case  of  exponential  tumor  growth  the  optimal  screening  schedules  with  a 
fixed  number  of  exams  have  a  trend  to  uniformity. 


Keywords  -  Screening,  Optimal  schedules,  Tumor  onset,  Tumor  size,  Carcinogenesis 
models 


1.  INTRODUCTION 

Because  of  the  significant  cancer  incidence  and  progress  of  tumor  detection  tech¬ 
nology,  cancer  surveillance  and  screening  are  becoming  increasingly  important  and 
costly  public  health  problems.  It  is  clear  that  appropriate  mathematical  methods 
are  indispensable  for  a  more  effective  management  of  the  caseload  through  designing 
optimal  surveillance  strategies.  Interest  in  exploring  this  avenue  has  quickened  in 
the  past  few  years  [1-17]. 

The  present  work  discusses  methodological  aspects  of  a  new  approach  to  opti¬ 
mization  of  cancer  screening  allowing  for  cancer  detection  at  the  earliest  stages  of 
tumor  development.  This  makes  the  chances  of  tumor  cure  more  favorable,  reducing 
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the  probability  of  tumor  recurrence.  The  problem  of  optimal  cancer  surveillance  is 
set  up  as  a  search  for  optimal  scheduling  of  screens  subject  to  certain  constraints  on 
the  number  and  timing  of  medical  exams.  Problems  of  a  similar  nature  have  already 
been  addressed  in  the  literature.  Yakovlev  and  Tsodikov  [1]  have  developed  methods 
for  constructing  optimal  surveillance  strategies  based  on  the  minimum  delay  time 
criterion,  given  that  the  total  number  of  examinations  is  fixed,  see  also  [2].  They  used 
dynamic  programming  methodology  to  solve  the  associated  optimization  problem. 
Their  results  show  that  this  approach  holds  much  promise  for  further  practical  use. 
As  one  example,  the  current  practice  for  the  breast  cancer  post-treatment  surveil¬ 
lance  at  Curie  Institute  (Paris,  France)  is  to  examine  the  patients  once  per  semester 
for  the  first  4  years,  once  per  year  for  the  next  6  years,  and  once  every  2  years  for 
the  remaining  period.  For  this  strategy,  the  estimated  false  negative  rate  appears  to 
be  equal  to  0.2  with  the  mean  delay  of  the  recurrence  detection  4.1  months.  Taking 
advantage  of  a  previously  proposed  parametric  model  of  tumor  recurrence  [3],  the 
authors  constructed  the  optimal  strategy  that  provides  a  33%  reduction  in  the  delay 
time,  with  the  tests  that  comprise  the  optimal  surveillance  schedule  tending  to  be 
more  frequent  when  the  hazard  rate  for  the  time  to  tumor  is  high.  However,  there 
are  two  weak  points  in  this  approach.  First,  the  probability  of  tumor  detection  is 
assumed  to  be  independent  of  the  process  of  tumor  regrowth.  Second,  estimation  of 
the  tumor  onset  time  distribution  is  feasible  only  if  a  sample  of  diagnostic  times  pro¬ 
duced  by  a  discrete  surveillance  program  with  known  false  negative  rate  is  available. 
The  same  applies  equally  to  pre-diagnosis  screening  programs. 

An  alternative  approach  to  the  problem  is  to  minimize  the  average  cost  of  surveil¬ 
lance  accounting  for  both  examination  costs  and  costs  of  late  detection  [1],  [4-15]. 
Since  the  two  cost  constituents  are  linked  in  the  optimization  procedure,  the  cost- 
utility  approach  makes  it  possible  to  search  for  both  the  optimal  number  of  exami¬ 
nations  and  their  sequence  in  time.  However,  the  costs  of  late  detection  are  usually 
very  difficult  to  evaluate.  For  yet  another  optimization  criterion  based  on  the  power 
of  a  statistical  test  for  mortality  rates,  the  reader  is  referred  to  [16], 

Focusing  our  effort  on  possible  medical  rather  than  economic  benefits,  we  propose 
to  explore  a  new  approach  to  the  problem  which  is  based  on  tumor  size  at  detection. 
Tumor  size  is  one  of  the  most  clinically  significant  characteristics  of  tumor  maturity 
that  determines  largely  the  probability  of  both  spontaneous  and  screening  based 
tumor  detection.  This  approach  makes  it  possible  to  utilize  data  on  tumor  size  at 
detection  as  an  additional  source  of  information  on  the  natural  history  of  the  disease; 
some  readily  available  epidemiologic  data  obtained  from  the  control  population  in 
the  absence  of  screening  appear  to  be  sufficient  for  estimation  purposes.  Another 
advantage  of  this  approach  is  that  it  offers  a  natural  way  for  incorporating  the  stage 
of  tumor  progression,  where  cancer  detection  normally  occurs,  into  stochastic  models 
of  carcinogenesis.  The  proposed  model  of  tumor  progression  accommodates  a  wide 
range  of  deterministic  and  stochastic  laws  of  tumor  growth. 

As  a  measure  of  the  effect  of  screening,  we  propose  to  use  the  difference  between 
the  expected  tumor  sizes  at  detection  with  and  without  screening,  which  coincides 
with  the  Kantorovich  distance  [18-21]  between  the  distributions  of  the  corresponding 
random  variables.  The  structure  of  this  distance  allows  for  characterizing  the  net 
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effect  of  screening,  as  compared  to  that  of  spontaneous  detection. 

Further  advancements  of  the  proposed  approach  to  constructing  optimal  sched¬ 
ules  of  cancer  screening  will  hopefully  give  answers  to  the  following  questions  of 
major  theoretical  and  practical  importance: 

1.  Is  the  optimal  efficiency  of  screening  high  enough  to  warrant  its  implementa¬ 
tion? 

2.  What  is  the  relation  between  the  optimal  screening  schedules  and  their  effi¬ 
ciencies  for  the  criteria  based  on  the  tumor  size  and  the  expected  time  delay? 

3.  What  are  cancer  specific  patterns  of  optimal  screening  schedules? 

4.  What  is  the  impact  of  hypothesized  laws  of  tumor  growth  on  the  optimal 
screening  efficiency  and  the  pattern  of  the  optimal  examination  schedules? 

5.  What  are  quantitative  characteristics  of  the  initiation,  promotion,  and  pro¬ 
gression  stages  for  specific  cancers? 

The  structure  of  the  present  paper  is  as  follows.  In  Section  2,  we  describe  some 
models  of  the  natural  history  of  cancer  (including  cancer  latency  and  growth),  screen¬ 
ing  schedules,  and  cancer  detection.  Here,  we  also  formulate  basic  assumptions  and 
introduce  mathematical  formalism.  An  explicit  formula  for  the  efficiency  functional 
is  derived  in  Section  3.  Sample  numerical  calculations  and  analysis  of  their  results 
are  addressed  in  Section  4. 


2.  BASIC  NOTIONS 
2.1.  Models  of  carcinogenesis 

In  describing  the  natural  history  of  cancer,  the  process  of  tumor  development  can  be 
broken  down  into  three  stages.  These  stages  are: 

•  formation  of  initiated  cells; 

•  promotion  of  initiated  cells  resulting  in  appearance  of  the  first  malignant  clono- 
genic  cell; 

•  subsequent  growth  and  progression  of  malignant  tumor. 

The  duration  of  each  stage  of  carcinogenesis  is  thought  of  as  a  random  variable 
(r.v.).  In  our  sample  calculations  presented  in  Section  4,  we  use  a  two-parameter 
gamma  family  to  specify  the  distribution  of  the  length  of  the  first  two  stages  of 
carcinogenesis.  However,  more  elaborate  mechanistic  models  of  carcinogenesis  are 
available  to  describe  the  time  to  the  event  of  malignant  transformation.  We  provide 
two  examples  of  such  models. 

The  most  widely  accepted  model  of  tumor  latency  is  commonly  referred  to  as 
the  Moolgavkar-Venzon-Knudson  (MVK)  model  [22,  23].  This  Markovian  two-stage 
model  involves  four  parameters  that  refer  to  the  rates  of  initiation  of  target  stem  cells 
(that  is,  formation  of  primary  precancerous  lesions),  and  rates  of  division,  death  or 
differentiation,  and  malignant  transformation  of  initiated  cells.  It  was  first  pointed 
out  by  Heidenreich  [24]  and  subsequently  by  Hanin  and  Yakovlev  [25]  and  Heidenre- 
ich,  Luebeck  and  Moolgavkar  [26]  that  these  four  parameters  are  not  jointly  identifi¬ 
able  from  time-to-tumor  data.  In  the  case  of  constant  parameters,  all  triples  of  their 
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identifiable  combinations  were  described  at  length  in  [25].  In  the  latter  case,  the 
MVK  model  leads  to  the  following  explicit  formula  for  the  distribution  of  the  total 
duration  T  of  the  first  two  stages,  that  is,  of  the  time  from  the  birth  of  an  individual 
to  the  tumor  onset  [27,  28]: 


FT(t)  :=  Pr(T  >  t)  = 


(< a  +  b)eat 


p 


b  +  ae(a+6)t  ’ 


t  >  0. 


(1) 


Here  a,  b,  p  >  0  are  identifiable  parameters  of  the  model,  Ft  :=  1  —  Ft  is  the 
survivor  function  of  the  r.v.  T,  and  Ft  is  the  cumulative  distribution  function  (c.d.f.) 
of  the  r.v.  T. 

Another  model  of  carcinogenesis  was  proposed  by  Yakovlev  and  Polig  in  [29]. 
According  to  this  model,  the  hazard  function  <f>  of  the  time  T  of  tumor  latency,  which 
is  related  to  the  survivor  function  by 

FT(t)  =  e-ti*{s)ds,  t>  0,  (2) 


is  of  the  form 

<f>(s)  =  9 ie"02 So  h^du  f  h(u)f{s  -  u)du,  s  >  0,  (3) 

where  h  is  a  given  time-dependent  rate  of  external  exposure,  /  is  the  probability  den¬ 
sity  function  (p.d.f.)  of  the  tumor  promotion  time,  and  9\,  62  are  positive  constants. 
The  key  feature  of  the  Yakovlev-Polig  model  is  that  it  allows  for  the  process  of  cell 
death  to  compete  with  the  process  of  tumor  promotion.  Two  particular  cases  of  the 
model  referring  to  spontaneous  and  induced  carcinogenesis  were  employed  in  [30] 
and  [31]  to  study  the  distribution  of  tumor  size  under  a  threshold  type  mechanism 
of  tumor  detection.  Recently,  Hanin  and  Boucher  [32]  found  conditions  under  which 
the  parameters  /,  0X,  02  of  the  model  given  by  (3)  are  identifiable  from  time-to-tumor 
observations.  Specifically,  a  general  necessary  condition  for  identifiability  of  model 
(3)  is  given  by  the  following  theorem. 

Theorem  1.  Suppose  that  the  function  h  satisfies  f£°  h(t)dt  <  00  and  that,  for  some 
C  >  0,  h(t)  =  0  for  t  >  C.  If  the  model  is  identifiable  in  a  family  T  then 

F{C)  >  0  for  all  F  €  T. 

Definition.  A  family  T  of  absolutely  continuous  probability  distributions  on  R+  is 
said  to  be  graduated  if  for  every  two  distinct  p.d.f ’s  /,  /  6f  and  for  every  constant 
A  >  0,  there  is  a  number  r  >  0  (which  may  depend  on  /,  /,  and  A)  such  that  either 
Af(t)  >  f(t)  for  all  t  >  r  or  Af(t)  <  f(t)  for  all  t>r. 

The  following  result  generalizes  Theorem  1  in  the  case  of  graduated  families. 
Theorem  2.  Suppose  that  h  is  bounded,  supported  on  [0,  C\  for  some  C  >  0,  and 
positive  almost  everywhere  on  [0,(7].  Then  the  model  is  identifiable  in  a  graduated 
family  T  if  and  only  if  F(C )  >  0  for  all  F  €  T . 
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2.2.  Tumor  growth 

The  following  general  functional  form  is  assumed  for  the  tumor  size  (the  number  of 
cells  in  a  tumor)  S  : 

S(w )  =  fe(w),  (4) 

where  w  is  the  time  from  the  moment  of  the  onset  of  cancer,  and  9  is  a  parameter 
which  may  be  scalar  or  vector,  deterministic  or  random.  It  is  assumed  that,  for  every 
6,  fo  is  a  strictly  monotonically  increasing  absolutely  continuous  function  such  that 
fd( 0)  =  1.  For  a  given  9,  denote  by  g&  the  inverse  function  for  fo,  and  set 

nw 

$e(w)  :=  /  fe(u)du. 

Jo 

Specific  laws  of  tumor  growth  of  primary  interest  are: 

(1)  Deterministic  exponential  growth;  in  this  case,  S(w)  =  eXw,  where  A  >  0  is  a 
constant  growth  rate,  see  [33]  for  substantiation; 

(2)  Exponential  growth  with  A  thought  of  as  a  gamma  distributed  r.  v.  [34]; 

(3)  The  Gompertz  law: 

S(w)  = 

with  constant  parameters  A,  B  >  0. 

2.3.  Screening  schedules 

The  sequence  of  moments  of  time  assigned  for  medical  exams  for  a  specific  cancer 
and  counted  from  the  birth  of  a  patient  will  be  called  a  screening  schedule.  Let  T  be 
the  set  of  all  possible  screening  schedules  r  =  {t\  <  r2  <  ...  <  rn}.  The  set  T  may 
be  subject  to  (some  of)  the  following  restrictions: 

(a)  n  <  no  i  where  no  is  an  upper  bound  for  the  number  of  exams; 

(b)  t i  >  rn  and  rn  <  M,  where  m  and  M  are  the  earliest  and  the  latest  times  for 
the  first  and  the  last  exams,  respectively; 

(c)  ri+ 1  -  n  >  h  >  0  for  all  i  =  1,2,  ...,n  -  1.  This  condition  suggests  a  lower 
bound  h  for  the  minimal  duration  between  any  two  successive  exams. 

Other  restrictions  on  the  moments  of  exams  can  also  be  accommodated.  In  the 
language  of  control  theory,  the  set  T  is  referred  to  as  the  set  of  admissible  schedules. 

2.4.  Tumor  detection 

We  distinguish  between  spontaneous  and  screening  based  tumor  detections.  The 
first  occurs  in  the  absence  of  or  concurrently  with  screening  and  is  thought  of  as  a 
continuous  process.  In  contrast  to  this,  screening  based  detection  is  an  instantaneous 
event  that  may  occur  only  at  the  moments  of  the  prescribed  medical  exams  and  is 
therefore  a  discrete  process.  When  both  types  of  detection  are  present,  they  can  be 
viewed  as  competing  risks. 

Numerous  attempts  have  been  made  to  relate  the  probability  of  detecting  a  tumor 
to  its  size  [33-37].  Following  Brown  et  al.  [37],  we  assume  that  the  rate  r0  of 
spontaneous  tumor  detection  is  proportional  to  the  current  tumor  size: 

r0  =  a0S,  (5) 
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where  ao  is  a  positive  constant. 

Let  r.v.’s  W0  and  W\  denote  the  times  of  spontaneous  and  screening  based  detec¬ 
tions,  counted  from  the  moment  of  cancer  onset,  respectively.  Then  for  the  moment 
W  of  combined  detection,  when  both  detection  mechanisms  are  in  place,  we  have 
W  =  min(Wo,  Wi).  Denote  by 

N0  —  fe(W0)  and  N  =  fe{W)  (6) 

the  corresponding  tumor  sizes  at  spontaneous  and  combined  detection. 

Keeping  in  mind  relation  (2)  between  the  survivor  function  of  an  absolutely  con¬ 
tinuous  nonnegative  r.v.  and  its  hazard  rate,  we  derive  from  (5)  that,  in  the  case  of 
non-random  parameter  6, 

FWo(w)  =  So  r°(u)du  _  g— «o  J”  fo(u)du  _  e-ao<t>e(w)'  (7) 


Therefore, 

FNo(n )  =  PwMn))  = 

and  hence 

EN0  =  1  +  J™  FNo(n)dn  =  1  +  J™  e~a^e^n))dn  =  1  +  j™  e~a°^  f'e(u)du.  (8) 

If  0  is  a  r.v.  then  an  additional  integration  in  (8)  with  respect  to  the  distribution  of 
9  is  required. 

In  particular,  for  non-random  exponential  tumor  growth  with  rate  A,  we  have 


FWo(w)  =  e-^eXw-1\  w>0, 

(9) 

FNo(n)  =  n>  1, 

(10) 

tn 

II 

I—1 

+ 

(11) 

ol  o 


Equation  (10)  suggests  that  in  this  case  the  r.v.  N0  has  a  translated  exponential 
distribution  with  parameter  a0/X-  If  A  is  a  r.v.  which  is  gamma  distributed  with 
parameters  p,  u,  then  it  follows  from  (11)  that 

ENq  =  1  ■+ - • 

a0u 

We  now  specify  the  distribution  of  the  r.v.  W\.  Recall  that  W\  is  the  time  of 
screening  based  detection  (in  the  absence  of  spontaneous  detection)  counted  from  the 
moment  of  appearance  of  the  first  malignant  clonogenic  cell.  Indeed,  the  distribution 
of  Wi  depends  on  the  selected  screening  schedule  r  =  {ti  <  r2  <  ...  <  Tn}.  For  the 
sake  of  convenience,  set  r0  :=  0  and  rn+ 1  :=  oo.  It  suffices  to  define,  for  every  t  >  0, 
the  conditional  distribution  of  W\  given  that  T  =  t. 

Let  Ti  <  t  <  Ti+i ,  0  <  i  <  n.  For  0  <  i  <  n  —  1  and  i  +  1  <  k  <  n,  define 
the  probability  pt(k )  :=  Pr(Wi  =  rk  -  t\T  —  t )  of  tumor  detection  at  the  k- th 
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screen  given  the  cancer  onset  at  moment  t,  and  by  Pt(oo)  =  1  —  J2k=i+iPt(k)  the 
corresponding  conditional  probability  that  tumor  is  not  detected  by  screening. 

We  introduce  a  discrete  analogue  of  the  hazard  rate  for  the  screening  based  de¬ 
tection  by 

!H  =  £  rt(k)K->,  (12) 

k=i-\- 1 

where  6X  stands  for  the  Dirac  measure  at  x,  and  the  sum  over  the  empty  set  of  indices 
is  set,  as  usual,  to  be  zero.  By  definition,  the  discrete  measure  is  related  to  the 
conditional  survivor  function  of  W\  given  that  T  =  t  through  the  equation 

iW-<M=e- w  >0,  (13) 

compare  with  (2).  It  follows  from  (12)  and  (13)  that 

n  n 

Y  Pt(j)+Pt(°°)  =  [J2Pt{j)+Pt( oo)]e~rt(fc) 

j=k+ 1  j—k 

or,  equivalently,  that 

1  -  Y  pt(j)  =  i1  -  Y  pt(j)}e~n{k) ■  (14) 

j=i+ 1  j=i+ 1 

For  k  =  i  +  1,  we  find  from  (14)  that 

1  —  pt(i  +  1)  =  e~rt(i+1\  (15) 

More  generally,  iterating  this  argument  we  obtain  that 

pt(k)  =  e~^=Y  rt(j)[l  -  e-rt(fc)],  i  +  1  <  k  <  n. 

Observe  that  this  holds  true  for  all  k  =  l,...,n,  if  we  set  Pt(k)  =  rt(k)  =  0  for 
1  <  k  <  i. 

Similar  to  (5),  we  are  assuming  that  the  discrete  rate  of  screening  based  detection 
is  proportional  to  the  current  tumor  size: 

rt(k)  =  aS(n-t),  i  +  1  <  k  <  n,  (16) 

with  some  constant  a  >  0.  Combining  (13),  (12)  and  (16)  with  (4)  we  find  that, 
given  any  t  such  that  r*  <t  <  Tj+i,  0  <  i  <  n  —  1, 

FwAT=t{w)  =  e-“^=<+1*^Tfc-t\  where  Tj  -  t  <  w  <  rj+1  -  t,  i  +  1  <  j  <  n. 

(17) 

Consider  the  case  of  one  exam  occurring  at  a  moment  r  with  the  detection  proba¬ 
bility  p  =  p(t,  t )  and  the  discrete  detection  rate  r  =  r(t,  r).  Then  by  (15),  1— p  =  e~r. 
If  the  probability  p  is  small  then  the  rate  r  is  approximately  equal  to  p.  In  partic¬ 
ular,  under  the  assumtion  (16),  the  probability  of  tumor  detection  is  approximately 
proportional  to  the  current  tumor  size:  p  ~  aS(r  —  t ).  Klein  and  Bartoszynski  [34] 
proceeded  in  their  study  of  breast  cancer  from  a  more  general  assumption  that  the 
probability  of  tumor  detection  is  proportional  to  some  power  of  the  tumor  size.  Their 
estimate  of  this  power  leads,  however,  to  a  value  which  is  very  close  to  1. 
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3.  FORMULA  FOR  THE  SCREENING  EFFICIENCY  FUNCTIONAL 


We  proceed  from  the  following  two  biologically  natural  assumptions. 

1.  The  r.v.’s  W0  and  T  are  independent. 

2.  For  every  t  >  0,  the  r.v.  ’s  W\  and  W0  are  conditionally  independent  given  that 
T  =  t. 

The  first  assumption  claims  that  the  moment  of  spontaneous  tumor  detection  mea¬ 
sured  from  the  appearance  of  the  first  malignant  clonogenic  cell  is  independent  of 
the  prior  duration  of  tumor  latency.  The  second  assumption  reflects  a  technological 
(or  instrumental)  nature  of  both  detection  processes.  It  states  that,  given  the  mo¬ 
ment  of  cancer  onset,  the  two  times  W0  and  W\,  at  which  competing  events  of  the 
spontaneous  and  screening  based  tumor  detection  may  occur,  are  independent.  This 
statement  immediately  follows  from  the  assumption  that  both  detection  processes 
are  completely  determined  by  the  current  tumor  size  as  a  deterministic  function  of 
time. 

For  an  admissible  screening  schedule  r  €  T,  we  define  the  efficiency  functional 
as  the  Kantorovich  distance  dK(N0,N;T )  (see  [18],  [20],  [21])  between  the  tumor 
sizes  N0  and  N  at  spontaneous  and  combined  detection.  This  quantity  serves  as  a 
clinically  natural  measure  of  the  gain  resulting  from  screening.  It  is  well  known  [19], 
[20]  that 

/CO  _  _ 

|  FNo(n)  -  FN(n)  \  dn.  (18) 

It  follows  from  (7),  inequality  W0  >  W,  and  monotonicity  of  the  function  f()  that  the 
r.v.  Nq  stochastically  dominates  the  r.v.  N  :  F^0  >  F^.  This  leads  to  the  following 
alternative  expression  for  the  efficiency  functional: 

/oo  _  r°°  _ 

FNo(n)dn-J^  FN(n)dn  =  EN0  -  EN,  (19) 

where  E  stands  for  the  expectation. 

Suppose  that  parameter  6  is  non-random.  We  set  n  =  fo(w)  and  condition  upon 
the  r.v.  T  in  (18)  to  obtain 


poo  _  _  , 

d(N,  N0 ;  r)  =  /  |  FWo(w)  -  Fw(w )  |  fg(w)dw 

Jo 

poo  poo  _  _  , 

=  Jo  Jo  '  Fwo<yW^  ~  Fw\T=t(w)  I  fe(w)dwdFT(t), 

where  Fw\T=t  is  the  conditional  survivor  function  of  the  r.v.  W  given  that  T  =  t. 
Since  W  =  mm(W0,  W ) ,  it  follows  from  our  assumptions  1  and  2  that 

Fw o  ~  Fw\r=t  =  Fwq  —  FWq  FWl \r=t  —  Fw0  F wy  \r=t  ■ 


Therefore, 


nco  _  t 

FWl\T=t(w)FWo(w)f0(w)dwdFT{t).  (20) 
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Observe  that  if  T  =  t,  where  rt  <  t  <  ri+ 1,  0  <  i  <  n,  then  the  only  possible  values  of 
ther.v.  W\  are  Tj+i— t,  ...,rn+i— t.  More  specifically,  Wx  —Tj-t,  i  +  1  <j  <n,ifthe 
j  — th  exam  detected  a  tumor,  and  Wx  =  rn+i  —  t  =  oo  if  the  tumor  was  not  detected 
in  the  course  of  screening.  Therefore,  if  t  >  rn  or  t*  <  t  <  rI+1)  0  <  i  <  n  —  1,  and 
0  <  w  <  Tj+i  —  t,  then  FWl\T=t(w)  =  0.  This  allows  us  to  rewrite  (20)  in  the  form 

d(N,N0;T)  =  Y/  E  /  FWllT=t(w)FWo(w)fd{w)dwdFT(t). 

i=0  Jri  j=i+ 1  Jt3  4 

We  now  recall  the  explicit  expression  (17)  derived  above  for  the  function  \r=t , 

and  denote 

Ge{x)  :=  /  FWo(w)fd(w)dw,  x  >  0, 

Jx 

to  obtain  finally 


d(N,  N( 


0, 


r)-tr 

i= o  Jr* 


E  I1 

j=i+l 


-  e-“SU.*<"-‘)][G,(Tj  -  t)  -  G»(Ti+1  -  tJldFrW 


=  l)/n+1  y  t)dFT(t).  (21) 

i=0  ^Ti  j=i+ 1 

In  the  case  when  parameter  0  is  random,  the  right-hand  side  of  (21)  should  be 
integrated  additionally  with  respect  to  the  distribution  of  9. 

If,  in  particular,  fe(w)  =  eXw  with  a  constant  rate  A,  then  invoking  (9)  we  find 
easily  that 

Ge(x)  =  —  x  >  0. 

a0 

In  this  case  the  efficiency  functional  (21)  takes  on  the  form 


A  ^ 

£/_ 

j=*+i 

Observe  also  that  (19)  implies 


d(N,N0]r)  =  —  Yl  r  E 

a0  i=o  Jri  j=i+l 

(22) 


EN  =  EN0-d(N,N0-,T). 


This  allows  for  an  explicit  calculation  of  the  expected  tumor  size  at  combined  detec¬ 
tion  on  the  basis  of  formulas  (8)  and  (21). 

The  problem 

d(N,  N0]  t )  — »  max,  t  G  T,  (23) 

can  be  solved  by  exhaustive  search  with  some  simplification  arising  from  the  special 
form  of  the  dependence  of  the  functional  (21)  on  r.  A  question  of  practical  impor¬ 
tance  is  what  are  the  values  of  the  number  n  of  exams  for  which  the  problem  (23) 
is  computationally  feasible.  We  will  conclude  this  paper,  which  deals  primarily  with 
methodological  and  mathematical  aspects  of  the  problem  of  optimization  of  cancer 
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surveillance,  with  some  sample  calculations  with  prescribed  values  of  model  param¬ 
eters. 


4.  NUMERICAL  EXPERIMENTS 

It  was  assumed  that  the  time  T  to  tumor  onset  is  gamma  distributed  with  the  mean 
fj,  —  50  years  and  the  standard  deviation  a  =  20  years.  The  graph  of  the  c.d.f.  Ft  is 
shown  in  Fig.  1.  The  law  of  tumor  growth  was  taken  to  be  deterministic  exponential 
with  the  rate  A  =  1.6  years-1,  which  corresponds  to  the  tumor  size  doubling  time  of 
approximately  5.2  months.  The  rate  of  spontaneous  tumor  detection  was  assumed 
to  be  a0  =  0.03.  The  graph  of  the  survivor  function  FWo  given  by  equation  (9) 
is  presented  in  Fig.  2.  The  effect  of  one  exam  occurring  after  tumor  onset  with 
the  screening  based  tumor  detection  rate  a  =  0.1  is  shown  in  Fig.  3  featuring  the 
survivor  function  of  the  time  W  to  combined  tumor  detection. 

The  search  for  optimal  screening  schedules  and  optimal  screening  efficiencies  was 
conducted  for  a  fixed  number  n  of  screens  with  no  restriction  on  the  moments  of 
exams  and  for  various  values  of  a.  The  method  of  optimization  was  the  exhaustive 
search  with  the  step  0.25  years.  Parameter  values  fx  —  50  years,  A  =  1.6  years-1,  and 
a0  =  0.03  were  fixed  throughout  the  calculations.  For  n  =  10,  plots  of  the  rescaled 
optimal  screening  efficiency  d  with  a  =  20  years  versus  a  and,  for  a  =  0.1,  versus  a 
are  shown  in  Fig.  4  and  5,  respectively.  As  it  could  be  expected,  d  increases  with 
increasing  a  and  decreases  with  increasing  a. 

The  results  of  our  search  for  optimal  screening  schedules  with  n  =  10,  20  and 
with  several  values  of  a  and  a  are  given  in  Table  1.  For  the  reader’s  convenience, 
screening  schedules  are  represented  by  the  intervals  A;  :=  Tj  —  T;_i,  i  =  l,...,n, 
between  two  successive  exams.  For  all  cases  explored,  optimal  screening  schedules 
are  uniform  or  very  close  to  such. 

As  a  test  for  optimality  of  a  screening  schedule,  profiles  of  the  efficiency  functional 
(22),  with  n  —  1  moments  of  exams  fixed  at  the  optimal  values  and  the  remaining 
one  varying  between  the  two  fixed  neighboring  moments  of  exams,  were  computed. 
For  n  =  20  and  a  =  0.1,  these  profiles  are  given  in  Fig.  6.  For  the  moment  T\, 
a  clear  cut  maximum  was  observed  (see  Fig.  6d),  while  for  t2 o  the  maximum  is 
more  flat  (see  Fig.  6c).  All  intermediate  moments  of  exams  r2,  ...7i9  demonstrated  a 
well- pronounced  parabolic  maximum  (see  Fig.  6a,  b). 
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Table  1. 


Optimal  screening  schedules 


n=20 


a 

a 

Ai=x, 

a2 

a3 

a4 

A5 

a6 

A  7 

As 

a9 

Aio 

20 

0.1 

27.75 

2.25 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.25 

A]} 

Al2 

An 

A14 

Aj5 

Ai6 

A|7 

Al8 

A|9 

o 

CS 

<1 

2.25 

2.25 

2.25 

2.25 

2.25 

2.25 

2.25 

2.25 

2.25 

2.25 

n=10 


a 

a 

A|— T] 

a2 

a3 

A4 

As 

A6 

A? 

As 

a9 

Aio 

20 

0.1 

35.00 

2.50 

2.50 

2.50 

2.25 

2.50 

2.50 

2.50 

2.50 

2.50 

20 

0.3 

35.00 

2.50 

2.50 

2.50 

2.50 

2.50 

2.50 

2.50 

2.50 

2.50 

20 

0.5 

34.50 

2.75 

2.75 

2.50 

2.50 

2.50 

2.50 

2.50 

2.75 

2.75 

10 

0.1 

40.75 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.25 

30 

0.1 

27.25 

2.50 

2.75 

2.50 

2.50 

2.50 

2.50 

2.50 

2.50 

2.50 
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